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Abstract We consider the potential-model approach for ob¬ 
taining the spectrum of charmonium and bottomonium, re¬ 
placing the usual gluon propagator by one obtained from lat¬ 
tice simulations. The resulting spectra are compared to the 
corresponding ones in the Cornell-potential case. We also 
estimate the interquark distance in both cases. 

1 Introduction 

A reliable description of heavy quarkonia states is of 
great interest for our understanding of nonperturbative as¬ 
pects of QCD [ 1 ] and is expected to be important in guiding 
the search for physics beyond the standard model [2]. A for¬ 
tuitous advantage in the study of such states is that, due to 
the large mass of the heavy quarks, various approximations 
may be adopted. For example, an expansion in inverse pow¬ 
ers of the heavy-quark mass m is performed in potential non- 
relativistic QCD (pNRQCD) [3], and lattice simulations (es¬ 
pecially for bottomonium systems) are applied to effective 
actions obtained by an expansion in powers of the heavy- 
quark velocity v/c. Similarly, in the relativistic quark model 
with the quasipotential approach, radiative corrections may 
be included and treated perturbatively in the case of heavy 
quarkonia [4]. This possibility of exploring different scales 
of the problem separately is also helpful in methods more di¬ 
rectly based on QCD, such as studies of Dyson-Schwinger 
and Bethe-Salpeter equations [5]. 

An early but still successful approach to describe heavy 
quarkonia is given by nonrelativistic potential models, to 
which relativistic corrections may also be added [6]. * The 
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'Note that these corrections may be computed from lattice data for the 
Wilson loop [7, 8]. 


idea is to view confinement as an “a priori” property of 
QCD, modeling the interquark potential to incorporate some 
known features of the interaction at both ends of the energy 
scale. The simplest such model, the Cornell — or Coulomb- 
plus-linear — potential [9-11], is obtained by supplement¬ 
ing the high-energy (perturbative) part of the potential with 
an explicit confining term. The resulting expression is a sum 
of two terms: the first one comes from the quark-antiquark 
interaction in the one-gluon-exchange (OGE) approxima¬ 
tion and the second one is a linearly rising potential. We 
have 

( 1 ) 

3 r 

where is the strong coupling constant and a is the string 
tension. The first term may be associated with scattering 
of the quark-antiquark pair inside the meson and is analo¬ 
gous to the Coulomb potential in the QED case. The second 
term corresponds to linear confinement as observed from 
the strong-coupling expansion of the Wilson loop in lattice 
gauge theory with static quarks. 

In practice, the static interquark potential may be defined 
conveniently in terms of the Wilson loop, or it may be ob¬ 
tained (perturbatively) by taking the nonrelativistic limit in 
the Bethe-Salpeter equation describing the bound state of 
two heavy fermions. This yields a Schrodinger equation, to 
which a linear term is added a posteriori. The numerical 
procedure for obtaining the mass spectrum for the Cornell 
potential, as well as for other commonly used potentials is 
reviewed in detail in [12]. 

The Cornell potential provides a spin-independent de¬ 
scription of the interquark potential for heavy quarks, with 
parameters determined by fitting a few known states (see e.g. 
[13]) or by comparison with lattice simulations. Eor a recent 
determination of these parameters, see Ref. [14]. It would 
be interesting, nevertheless, to have a better insight about 
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confinement as an emergent property of the interquark in¬ 
teraction induced by the gluon propagator, rather than as a 
built-in feature. 

Of course, the gluon propagator in QCD is very differ¬ 
ent from the perturbative one at the hadronic scale, and it 
should contain full information about confinement. In order 
to use this nonperturbative information we propose to sub¬ 
stitute the free gluon propagator in the OGE term of the po¬ 
tential, as described above, by a fully nonperturbative one, 
obtained from lattice simulations. We want to check if this 
replacement leads to an improved description of the spec¬ 
tra, possibly without the need to include the linearly rising 
term explicitly. We use the data generated in studies of the 
SU(2) gluon propagator in Landau gauge on very large lat¬ 
tices (up to 128'*), reported in [15, 16]. More precisely, we 
use directly the fit obtained in Ref. [17]. We note that our 
aim is to gain a qualitative understanding of the interplay 
between perturbative and nonperturbative features of the in¬ 
terquark potential. Our approach is similar in spirit to the 
one in Refs. [18, 19], but our conclusions are different. 

Preliminary versions of our study have been presented 
in Refs. [20] and [21]. 

We organize this paper in the following way. In Section 
2 we review the procedure for obtaining the Coulomb po¬ 
tential in QED as the nonrelativistic limit of scattering 
and how this is adapted to heavy quarks. In particular, we in¬ 
troduce the lattice propagator in the OGE term. In Section 3 
we describe our method for obtaining the mass spectra with 
the desired potential. Our results are presented in Section 4 
and our conclusions in Section 5. 


2 Potential from Lattice Propagator 

Let us hrst review how the Coulomb potential is ob¬ 
tained in the nonrelativistic limit of QED from the applica¬ 
tion of Eeynman rules to the electron-positron system. The 
scattering-matrix 5/,, from which the interaction potential 
may be obtained, is given by 

Sfi = {f\i) = dfi + i(2Kf 5(4) (Q - P) Tfi , (2) 

where Q and P correspond respectively to the final and ini¬ 
tial total momentum and 7y, is the scattering amplitude. The 
two tree-level Eeynman diagrams contributing to 7y, (see 
Eig. 1) correspond to the t and s channels, respectively com¬ 
ing from scattering with one photon exchange and to anni¬ 
hilation and creation of an pair. We get 

1 

“ pi? 

where 

^exch : ^1) 

X v(p 2 , 0 - 2 ) 7 ''(4) 



Fig. 1 Feynman diagrams corresponding to the two terms in the 
scattering amplitude. The left diagram corresponds to the t channel 
(photon exchange) and the right diagram to the i channel (pair annihi¬ 
lation). 


and 

Lnnihil — C v(p2) O2) T*'^(Pl i Oi) P^v(^) 

X m(^1,Ti)7^v(^2,T2) . (5) 


We follow the notation in [12, 22, 23]: p, denotes the mo¬ 
mentum of the incoming particles and qj of the outgoing 
ones. The particles’ initial and final spins are respectively 
(7, and T, . We represent the photon propagator by a function 
Pjiv (k) of the photon momentum k. 

We then make the nonrelativistic approximation, i.e. we 
impose the kinetic energy of the system to be much smaller 
than its rest energy (|p| m = £). The four-component state 
vectors become 



0 

0 


1 

0 

“(a 1/2) = 

; m(7’,-1/2) = 


Vo^ 




and 




0 

0 

vJ 




0 

-1 

Vo/ 


( 6 ) 


(7) 


In this approximation, Tfi can be written as 


Tfi = 


w 


■ (fexch T ^annihil) ■ 


( 8 ) 


To compute the exchange term, we adopt the Dirac represen¬ 
tation for the gamma matrices and the center-of-momentum 
frame, obtaining 


fexch = 5^^5aiXiPnv{k)5''^5a2t2 

= s 7*00 (^) 5(jj ti 5(J2 1 (9) 

with^ 


k = pi-q^= (0, k) . (10) 

^-We are using the metric g^v = diag(—1,1,1,1), which will be more 
convenient when we consider -Wick rotations later. 















3 


For the annihilation term, note that conservation of momen¬ 
tum at the vertices implies that 

k={2m,Q). ( 11 ) 


we obtain a Coulomb-like interquark potential 


y(r) 


4^ 
3 Anr 


4 

3 r ■ 


(19) 


For QED, the Feynman-gauge propagator is given by the 
expression Pjxv{k) = —gfiv/k^- As seen in Eqs. ( 6 ) and (7), 
the spinors are momentum-independent in the nonrelativis- 
tic approximation and, while fexch is proportional to 1 /k^, 
we see that fannihii will be proportional to 1 /4 ot^. Thus, we 
can neglect annihilation effects and the scattering amplitude 
is given by 


Tfi = 


1 ^2 
{2k)^ k 2 ■ 


( 12 ) 


The potential can then be obtained as an inverse Eourier 
transform, which leads to the Coulomb potential 


y(r) = —{27t)^J exp(—ik-r) Tfi{^)(P’k 




Notice that the above potential is non-confining. This 
is expected, since it results from a perturbative calculation, 
while confinement is a nonperturbative phenomenon. The 
addition of a linear term as described in Section 1 leads to 
the Cornell, or Coulomb-plus-linear, potential [9-11] 

y(r) = -f (7r, (20) 

3 r 

which describes surprisingly well the states of charmonium 
and bottomonium. 

As mentioned in Section 1, we substitute the free propa¬ 
gator by a fully nonperturbative one in the OGE term. More 
precisely, we use the propagator 


Eor QCD, we replace the photon by the gluon and the 
electron-positron pair by a quark-antiquark pair. The scat¬ 
tering amplitude will continue to be expressed as a sum of 
the two terms, now given by 

fexch = glu{q]_,Xl)c\jX“'f C]_ju{pi,ai) P‘^y{k) 

X v(/72, 0-2)4 C2,/v(^2,T2) (^4) 

and 

f a nnihii = - go v(p2, O2) 4 Cl, m(pi , (7i ) F"* (k) 

X m(?1,Ti)c4^A*7''c2,/v(^2,T2), (15) 

where C(i 2 ),(;,/) are three-component color vectors and A" 
are the Gell-Mann matrices. 

With respect to Eq. (3), the terms fgxch and fannihii will 
have multiplicative (Casimir) factors, coming from the sum 
over colors. This sum is obtained assuming that the incom¬ 
ing/outgoing quarks and antiquarks have equal probability 
of being in a given color state and imposing a color-diagonal 
gluon propagator. The factors are given respectively by 

1 saa a 

4 ^AVi,,c 4 AV 2 ,/ = -TrA“A« = ^ = 3 (16) 

and 

cl.X‘^cucljX“c 2 j = i(TrA")(TrA«) = 0. (17) 

Therefore, annihilation effects do not contribute, inde¬ 
pendently of the nonrelativistic approximation. If we now 
assume a free (i.e. tree-level) gluon propagator 


pah _ Sl^v 5 "* 

- ^2 ’ 


( 18 ) 



C{s-\-k'^) / kpkyX 

t'^ + u^k^ + k^ y k^ ) 




C = 0.784, s = 2.508 GeV4 

u = 0.768GeV, t = 0.720GeV4 


( 21 ) 


obtained from fits of lattice data for a pure SU (2) gauge the¬ 
ory in Landau gauge given in Ref. [17]. We note that the lat¬ 
tice data for propagators in the SU (2) and SU (3) case have 
essentially the same behavior apart from a global constant 
[24]. Also, the above parameters correspond to a value 1 /k^ 
at 2 GeV. Here we choose to normalize the propagator to 
1 jk^ at k —> 00 , i.e. we adopt C = 1. 

We now follow the same procedure as in the QED case. 
Erom Eq. (9) we notice that, in the nonrelativistic approx¬ 
imation, only the component Foo( 0 ,k) survives in the fexch 
term and thus the term kpky/k^ vanishes [see Eq. (10)]. 
Lastly, in order to convert the propagator in Eq. (21), which 
was obtained in Euclidean space, to Minkowski space, we 
undo the Wick rotation, taking 5pv ~gtiv- We obtain^ 



f^ + M^k^ -|- k^ 


( 22 ) 


This leads us to the following scattering amplitude 


^ 4 C(. + k4 

3 f^-f M^k^-f k^ 


(23) 


The potential is obtained, as was done in the QED case, 
as a Eourier transform of the scattering amplitude [see Eq. 

^Let us recall that the propagator is a gauge-dependent quantity. 
A gauge-independent potential obtained from the (Coulomb-gauge) 
propagator is discussed in [25]. 











4 


(13)]. We use spherical coordinates for k and set r = rt The 
angular integration is then trivial, resulting in^ 


y(r) 


4 sic 

3 2i{2nfr 


{s + k^) 

+ u^k? +k^ 


kdk. (24) 


The integral in Eq. (24) can be solved using residue calcu¬ 
lations. The four poles in the integrand are symmetrically 
distributed in the four quadrants of the complex plane. We 
index these poles in the following way 


km,n = (-l)'”!Vfexp 


(-iry 


OT, n = 0,1. 


(25) 


where 


9 


= arctan 



(26) 


Fig. 2 Comparison between the lattice-gluon-propagator potential 
TlGP the Coulomb-like potential (color factor included) in the 
charmonium case. 
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The associated contour integral is performed by considering 
its two terms separately; for the term with (respectively 
with we close the contour above (respectively below). 
The residues are given by 


Also, the nonrelativistic approximation removes any spin 
dependence from the interactions. This means that, in our 
description, states with different spin values will be degen¬ 
erate. 


Res 


(s + k^)ke^‘^'' 


1 {s + kl„) 

2 u^ + 2A:2 ^ 11 3 Method for Obtaining Quarkonium Masses 


The result is simplified by noticing that ki Q = —kQ Q, k^ i = 
—ki I and kw = k^^. The obtained potential, which we call 
the lattice-gluon-propagator potential Vlcp^ is given by 


VLGp{r) = 


C{s + kl^)e‘'^-o'- 

m2 -p 2^^ q 


(28) 


where = g^/AK [see Eq. (19)]. 


In order to use the above expression in our spectrum cal¬ 
culation, we use the four-loop formula available in Ref. [26, 
Section 9] as well as the values (also available in Ref. 

[26, Section 9]) to evaluate the strong coupling constant a, 
at the energy scale of the mass of the IS quarkonia states 
[respectively J/\j/ and r(lS) in the charmonium and bot- 
tomonium cases]. The resulting potential is compared (for 
the charmonium case) to the Coulomb-like potential from 
Eq. (19) in Eig. 2. We see that, although the two curves 
are clearly different, Vlgp is also non-confining. Thus, since 
(tree-level) perturbation theory was applied to obtain this 
potential, the property of confinement was lost, even though 
the used propagator was obtained nonperturbatively. 

We therefore add a linearly rising term or to the poten¬ 
tial in order to model confinement, as done for the Cornell- 
potential case. The resulting expression is the lattice-gluon- 
propagator-plus-linear potential 


VLGP+L{r) = VLGp{r) + Or. (29) 

^For the evaluation of this integral only, we will denote |k| =k. 


Consider a central (nonrelativistic) potential describing 
the interaction between two particles. Since we are dealing 
with a two-particle system, we can write the Hamiltonian 
in terms of relative coordinates and use separation of vari¬ 
ables in the resulting partial differential equation to isolate 
the angular part of the wave function, given by the spherical 
harmonics. Lastly, we perform the usual substitution of vari¬ 
ables in the radial wave function R{r) = /(r) jr to obtain the 
ordinary differential equation (ODE) for /(r) 


dr^ 


2ji 


E — V{r) — 2m 


2/4 r2 


fir) 


= 0 , 


(30) 


where ji is the reduced mass 
m 



(31) 


and m is the mass of the heavy (charm or bottom) quark. We 
use units such that c = h= 1. Notice as well the addition of 
the rest mass of the particles, which will allow us to compare 
the eigenvalue directly with the mass values given in Ref 
[26]. 

The above ODE has to be solved with proper boundary- 
value conditions. The first one is that /(O) = 0. This comes 
from the requirement that R{0) be non-singular. A second 
condition is that /(r —> oo) = 0 and comes from the fact that 
R{r) is normalized, i.e. 


r\R{r)\^r^dr = H \fir)\^dr = I. 
JO Jo 


( 32 ) 




























































5 


In the limit of large r, the potential is dominated by the 
linearly rising term and the ODE becomes 


numerical method to compute the eigenenergies for a differ¬ 
ent set of parameters, allowing comparison with our results.^ 


fir') 

^y-2^arf{r) =0. (33) 

The general solution of this equation is the linear com¬ 
bination of the Airy functions Ai{p) and Bi{p) [27], where 
p = r. However, the Airy function of the second 

kind Bi(p) diverges at large p and therefore it does not obey 
the boundary condition at infinity. For p > 0, the Airy func¬ 
tion of the first kind can be written as 

where Ky (x) is the modified Bessel function of the second 
kind. One can try the Ansatz /(p) = g(p) Ai(p) and use the 
property Ky (x) = V Ky (x) /x — Ky+i in the ODE in Eq. (30) 
to obtain a second-order ODE with coefficients in terms of 
Ai(p) and K 4 / 3 (x). However, by expressing these functions 
as a power series in p, one clearly sees that an analytic solu¬ 
tion would be challenging, even for the simpler case of the 
Cornell potential. We therefore seek a numerical solution of 
the problem. 

To this end, we use the so-called shooting method [28]. 
It consists in picking trial values in a discretized range for 
the eigenenergies, integrating the ODE for each of these val¬ 
ues to obtain the corresponding wave function and choosing 
the energies for which the wave function obeys the bound¬ 
ary conditions approximately. We use the backward second- 
order Runge-Kutta method to integrate the wave function, 
starting from a maximum value r^ax for the radial coordi¬ 
nate until the origin, in steps of lir (we adapt the method 
as presented in Ref. [28] by adopting a negative integration 
step). We choose r^ax sufficiently large so that we can use 
f{r,nax) = Ai{fmax) and f{rmax) = Ai’{rmax) as initial condi¬ 
tions. In practice, the wave function will not obey the bound¬ 
ary conditions exactly since the proposed energy is unlikely 
to be an exact eigenenergy. Nevertheless, we may count the 
number of nodes of the wave function: each time we observe 
an increase in the number of nodes when compared with 
the previously proposed energy, the desired eigenenergy will 
be between the two proposed values. We further refine our 
method by adapting the bisection method to search for the 
eigenenergy in this interval, thus allowing the use of a coarse 
grid without loss of precision. 

We test our algorithm with the parameters of the Cor¬ 
nell potential set to (7 = 1 GeV^, 2p = 1 GeV and Aas/lt = 
1. These parameters are the ones used in Ref. [29], which 
adopts a different approach for solving the problem. We find 
agreement with their values up to the 4th and in some cases 
even 5th decimal place. Similarly, Ref. [13] uses yet another 


We consider as free parameters in the ODE the string 
tension CJ and the mass of the quark. The reason for includ¬ 
ing the quark mass as a parameter is that quark masses are 
not observable directly and depend on the renormalization 
scheme. Therefore, we have the freedom to pick one that 
best describes the observed spectrum. 

To find the best values for these parameters, we set up a 
two-dimensional grid of values of m and a. Then, we com¬ 
pute the eigenenergies for each proposed set of parameters 
and select the one that best describes the observed spectrum. 
As a criterion for choosing the optimal parameters we con¬ 
sider the minimization of in the description of a few input 
values from experiment, i.e. we pick the set of parameters 
minimizing 

X"(parameters) = ^ ^ 

where CJ,- is the experimental error associated with the energy 
Ei.experimentai and the E/’s are the eigenenergies computed as 
described above. 

Notice that the numerical method described above can 
be applied to any central potential. Therefore, we can obtain 
the results for the Cornell and the Viqp+l potential using 
the same algorithm, with the same inputs, making it easy to 
compare the differences between the two cases. 


4 Results 


We apply our method to two closely nonrelativistic sys¬ 
tems: charmonium and bottomonium. Since the bottomo- 
nium is heavier in comparison with its kinetic energy, we 
expect that the proposed model will work better for it than 
for charmonium. 

As stated at the end of Section 2, this model does not 
see spin interactions. This results in states with high degen¬ 
eracy, in comparison with the experimental data. In order to 
perform a fit to tune the parameters of our model, we need to 
average over these states with different spin. This is done in 
Ref. [30] by using the degeneracy of each state as a weight, 
i.e. the spin-averaged mass of the states with principal quan¬ 
tum number n and in the A-wave state {X = S,P,D...) is 
given by 


{M{nX)) 


Lfii mi{nX)gi 



(36) 


^More precisely, with respect to our Tables 3, 4 and 5 in Section 4, 
we find agreement up to the 3rd decimal places. We must consider 
that, in this comparison, our parameters are close to but not identical 
to the ones used in Ref. [13], which might explain the slightly worse 
agreement than in the comparison with [29]. 
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where mi{nX) is the mass of each of the Ni states with the 
same / and gi is the degeneracy of the state. The uncer¬ 
tainty associated with the above average may be estimated 
by propagation of errors, where the error of each mass is due 
to the width of the resonance peak®. This leads to very small 
errors, and we find that this is not a good way to estimate our 
uncertainties, as discussed below. We refer to this averaging 
procedure as “Mass Input 1”. 

Instead, we choose to average the spins by imagining 
that, if the experiments were not very precise, we would not 
see several narrow nondegenerate states, but broad degen¬ 
erate ones, i.e. a low precision experiment would see the 
peaks merged. We thus consider as input for a state the mid¬ 
point between the state with lowest energy and the one with 
highest energy. The error is estimated as half of the distance 
between these two states. We refer to this method as “Mass 
Input 2”. 

In Tables 1 and 2 we show the data extracted from Ref. 
[26] for charmonium and bottomonium, combining the dif¬ 
ferent spin states using the two methods described above 
(i.e. “Mass Input 1” and “Mass Input 2”). We choose to in¬ 
clude just the states present in the meson summary table of 
Ref. [26] that are regarded as established particles. We omit 
charged states from our table, since quarkonia states must 
be neutral. We remark that, for our fit, we use only the states 
IS, IP and 2S of charmonium and bottomonium. 

We implement the algorithm described in Section 3 us¬ 
ing the following parameters: for the charmonium case, 
the quark mass ranges from l.OGeV to 2.0GeV in steps 
of 0.01 GeV, the string-tension parameter a ranges from 
O.lGeV^ to 0.5 GeV^ in steps of 0.01 GeV^ and we use^ 
Us = 0.2663. The wave function is integrated from a max¬ 
imum distance of 20.0GeV^* until the origin in steps of 
5.0 X lO^^GeV^* and the eigenenergies are searched in the 
range from 2.0GeV to 6.0GeV in steps of 0.04GeV. For the 
bottomonium case, the quark mass ranges from 4.15 GeV 
to 4.85 GeV in steps of 0.01 GeV, we use the same range as 
above for CJ and = 0.1843. The parameters of integration 
for the wave function are the same as for the charmonium 
case and the eigenenergies are searched in the range from 
8.5GeV to 12.5GeV in steps of 0.04GeV. 

We first made fits for the charmonium and bottomonium 
spectra independently, using “Mass Input 1”. This yields 
a very large value of for both systems (varying from 
1.9 X 10^/dof — when we fit bottomonium with Vlgp+l — 
to 23 X lO^/dof — in the case where we fit charmonium 
using Cornell potential). We thus conclude that the errors 

®We recall that bound states are detected by plotting a histogram of 
number of particles (cross-section) detected in a collision versus the 
energy of the collision. When a resonance is found, it is associated to a 
bound state. 

^We actually compute and use ctj with 20 decimal places to ensure 
precision, but we are aware that the error of the computation must be 
far greater. 


Fig. 3 Comparison of the Cornell Potential with Wcf+l- Fur the 
value of the strong coupling constant, we choose the one used in the 
(constrained-fit) description of the charmonium spectrum. 



shown in the column “Mass Input 1” of Tables 1 and 2 are 
inconsistent with the overall precision of our calculation, 
considering e.g. the effects of the nonrelativistic approxima¬ 
tion employed. Instead, using “Mass Input 2” gives accept¬ 
able values for in the charmonium case, but large val¬ 
ues for bottomonium (although orders of magnitude smaller 
than using “Mass Input 1”). The situation improves if an 
unconfirmed state of bottomonium is included [namely, the 
77^(15)]. Let us note that, in all cases, using Vlgp+l yields 
similar results to the Cornell-potential case, with eigenstates 
slightly closer to the experimental values and with smaller 
errors. 

The data obtained in the independent fits of charmo¬ 
nium and bottomonium spectra are then used to make a con¬ 
strained fit, i.e. one with a common value for the string ten¬ 
sion a of the two systems^. Notice now that we will have 
three free parameters (m^, and a) and six states as in¬ 
puts in the fit (the states IS, 2S and IP of charmonium and 
bottomonium), resulting in three degrees of freedom. The 
results of this fit using the Cornell potential and Vlgp+l are 
shown in Tables 3 and 4. 

In order to establish a confidence level for our param¬ 
eters, we use the method described in detail in Ref. [28], 
which consists in determining the region in parameter space 
for which x^ increases by less than one unit with respect to 
its minimum value, for each of the parameters. The confi¬ 
dence level obtained for each parameter is shown in Table 
5. This also allows us to establish confidence levels for the 
eigenenergies. In cases for which the obtained confidence 
level is asymmetric, we adopt the larger value as the error. 
We then draw three random numbers, each one following a 
Gaussian distribution centered at the value of the optimal pa¬ 
rameter found and with standard deviation equal to the sym¬ 
metrized error. We use these numbers as parameters to com- 

*We consider as input to this constrained fit just our second input 
method, i.e. column “Mass Input 2” of Tables 1 and 2, including the 
unconfirmed ?7(,(15) state. 
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Table 3 Results for the charmonium eigenstates using Vicp+l and the 
Cornell potential in a constrained fit (see text). 


Chai'monium Spectrum 

Vlgp+l 

State 

Predicted mass 

Deviation from averaged 


(GeV) 

spin state (GeV) 

IS 

2.96(5) 

-0.10 

2S 

3.69(6) 

0.01 

3S 

4.22(7) 

- 

IP 

3.46(6) 

-0.07 

2P 

4.02(6) 

0.09 

3P 

4.50(7) 

- 

ID 

3.81(6) 

- 

2D 

4.31(6) 

- 

Cornell Potential 

State 

Predicted mass 

Deviation from averaged 


(GeV) 

spin state (GeV) 

IS 

2.93(14) 

-0.14 

2S 

3.69(17) 

0.01 

3S 

4.28(20) 

- 

IP 

3.42(16) 

-0.11 

2P 

4.04(18) 

0.12 

3P 

4.58(21) 

- 

ID 

3.80(17) 

- 

2D 

4.36(20) 

- 


pute the spectrum, repeating the process N = 1000 times. 
The obtained values are presented in parentheses in Tables 
3 and 4. The resulting potentials are shown for the charmo¬ 
nium case in Fig. 3. A visual representation of the spectrum 
is provided in Figs. 4 and 5. We also compare these results 
with the ones in Ref [13], obtaining agreement with their 
eigenenergies to within ~ 10^^. 

For both the charmonium and the bottomonium spec¬ 
tra, we obtain smaller errors and better agreement with the 
spin-averaged experimental values in the Vicp+l case than 
in the Cornell-potential one (see Tables 3 and 4). This is 
especially true for the charmonium spectrum. In particular, 
the charmonium spectrum obtained using the Cornell poten¬ 
tial is not in agreement with experiment. For bottomonium, 
instead, the comparison shows agreement. This is not sur¬ 
prising, since bottomonium is a more closely nonrelativistic 
system. Of course, it would be interesting to check if the 
inclusion of relativistic corrections would allow a better de¬ 
scription of the charmonium spectrum and higher accuracy 
in the bottomonium case. We point out that the errors in the 
Vlgp+l case are of the same order of magnitude as in the 
spin-averaged experimental values for most states. 

An advantage of our approach is that we have direct ac¬ 
cess to the radial wave function f{r). We plot, as an exam¬ 
ple, the wave functions^ for the IS state for both potentials 
in the charmonium and bottomonium cases in Fig. 6. We can 

^The wave functions obtained using our code are not normalized. We 
interpolate the data and normalize f{r) a posteriori. 


Table 4 Results for the bottomonium eigenstates using Vlgp+l and the 
Cornell potential in a constrained fit (see text). 


Bottomonium Spectrum 

Vlgp+l 

State 

Predicted mass 

Deviation from averaged 


(GeV) 

spin state (GeV) 

IS 

9.47(4) 

0.04 

2S 

10.00(4) 

-0.01 

3S 

10.37(5) 

0.01 

4S 

10.67(6) 

0.10 

5S 

10.95(7) 

- 

IP 

9.86(4) 

-0.03 

2P 

10.24(5) 

-0.01 

3P 

10.56(5) 

0.03 

4P 

10.84(6) 

- 

ID 

10.11(5) 

-0.05 

2D 

10.44(5) 

- 

3D 

10.73(6) 

- 

Cornell Potential 

State 

Predicted mass 

Deviation from averaged 


(GeV) 

spin state (GeV) 

IS 

9.49(8) 

0.06 

2S 

10.00(10) 

-0.01 

3S 

10.39(12) 

0.03 

4S 

10.72(14) 

0.14 

5S 

11.01(16) 

- 

IP 

9.84(10) 

-0.04 

2P 

10.25(11) 

0.00 

3P 

10.59(13) 

0.05 

4P 

10.89(15) 

- 

ID 

10.10(11) 

-0.06 

2D 

10.45(12) 

- 

3D 

10.77(14) 

- 


Table 5 Quark masses and string tension obtained from our fit. These 
parameters are used to obtain the spectrum in Tables 3 and 4. 


Vlgp+l 

Cornell Potential 

Quark Mass 
in Ref. [26] 

mc= 1.16(3) GeV 

mc= l.ll+;;;j;^GeV 

Me = 1.275(25)GeV 

m* =4.617^:1]? GeV 

nii = 4.58i[{;}]]GeV 

m(,(MS) =4.18(3)GeV 

(7 = 0.23(1) GeV 

fr = 0.26;“[GeV 

mj(lS) =4.66(3) GeV 

= 6.20 

21^ = 12.13 



see that the similarity between the two potentials (see Fig. 3) 
and the obtained spectra is present for the wave functions as 
well. Also, note that the wave function is more extended for 
the charmonium states. 

This direct access to the the wave function can be of in¬ 
terest in other applications, such as effective field theories, 
for which one needs information on the typical distance be¬ 
tween the quarks [1]. We estimate this quantity by comput¬ 
ing 

roo 

d = 2{r) = 2 rfirfdr. (37) 

Jo 




























Fig. 4 Experimental mass spectrum for charmonium, together with the spin averages used as input in our calculations. We also show our results 
in the Vlgp+l and Comeli-potential cases. The fit shown is that of the constrained case and considering as input the states IS, 2S and IP of both 
systems. 
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Fig. 5 Experimental mass spectrum for bottomonium, together with the spin averages used as input in our calculations. We also show our results 
in the Vlcf+l and Cornell-potential cases. The fit shown is that of the constrained case and considering as input the states IS, 2S and IP of both 
systems. 
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Fig. 6 Comparison of the wave function f[r) of the IS state for bot- 
tomonium and charmonium using both potentials. 



Table 6 Typical interquark distances for charmonium. Errors are ex¬ 
pected to be negligible. 


Charmonium 

State 

Vlgp+l 

Cornell Potential 


distance (fm) 

distance (fm) 

IS 

0.80 

0.83 

2S 

1.59 

1.56 

3S 

2.20 

2.14 

IP 

1.29 

1.28 

2P 

1.94 

1.90 

3P 

2.50 

2.43 


Table 7 Typical interquark distances for bottomonium. Errors are ex¬ 
pected to be negligible. 


Bottomonium 


State 

Vlgp+l 
distance (fm) 

Cornell Potential 
distance (fm) 

IS 

0.45 

0.47 

2S 

0.94 

0.94 

3S 

1.35 

1.31 

IP 

0.76 

0.77 

2P 

1.18 

1.16 

3P 

1.55 

1.49 


We note that the factor 2 in Eq. (37) comes from the fact 
that r corresponds to the distance between one of the quarks 
and the center-of-mass of the system. Some of these typical 
distances are presented in Tables 6 and 7. 

Finally, we could also estimate decay widths, which are 
proportional to |/?(0)p. Notice, however, that this calcula¬ 
tion would require a more strict control of the numerical 
integration in the region near the origin, since the function 
R{r) = f{r) jr typically shows a divergence for r —>■ 0. 


5 Conclusions 

We briefly reviewed the potential-model approach for 
determining the spectrum of quarkonia and discussed the 
simplest such approach, the Cornell potential. We then mod¬ 
ified the procedure by replacing the free gluon propagator 
with one obtained using lattice simulations. The resulting 
Vlcp potential is different from the Coulomb-like potential, 
but is still non-confining. Inspection of Fig. 2 shows that, 
up to the hadronic scale, the potential rises above zero, with 
a trend to rise further. This is no longer true for larger val¬ 
ues of r, for which the potential is damped. In fact, in order 
to obtain a confining (linear) potential, the gluon propaga¬ 
tor should show a strong divergence, of 1 /k^, in the infrared 
limit, as proven in Ref. [31]. Also, an oscillating behavior 
— due to the complex poles of the lattice propagator [17] — 
is observed. 

We therefore add a linear term to Vi^qp, obtaining the 
Vlgp+l potential in Eq. (29). We solve the associated Schro- 
dinger equation numerically and compare our results with 
the spin-averaged spectrum. This is done for the Cornell po¬ 
tential and for Vlgp+l, both for charmonium and for bot¬ 
tomonium states. The spectra computed using Vlgp+l show 
a slight improvement over the ones obtained using the Cor¬ 
nell potential, but no qualitative differences are observed. In 
particular, the resulting potentials are rather similar, as seen 
in Fig. 3. A more accurate description of the spectra could 
be achieved by introducing relativistic corrections, as was 
done in Refs. [6, 30]. 

We were also able to obtain the wavefunction for all the 
states mentioned, which allows us to estimate the interquark 
distance of the considered states. Let us note that the wave 
functions are remarkably similar for the Cornell potential 
and Vlgp+l (see Fig. 6), even though the potentials are not 
identical (see Fig. 3). This might suggest that the wave func¬ 
tion is somewhat insensitive to details of the potential. In 
fact, a visual comparison between our wavefunctions and 
the one presented in [32, Fig. 5] (corresponding to a differ¬ 
ent parametrization of the Cornell potential) shows that they 
are also similar. 

Let us mention that a study using a similar method was 
carried out in Refs. [18, 33] to propose a potential for heavy- 
quarkonium states. In that case, the gluon propagator was 
taken from a study of Schwinger-Dyson equations [34]. This 
propagator is in qualitative agreement with the lattice re¬ 
sults we use. The main difference with respect to our study 
is that these authors do not include the linear term in the 
potential, but consider an additive contribution to the one- 
gluon-exchange potential, in such a way that the zero of the 
proposed potential coincides with the Cornell one. This cor¬ 
responds to fixing the self-energy of the static source [35], 
which appears in higher-loop calculations. A constant term 
in the interquark potential can also be interpreted as being 
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related to the infrared divergence of the Fourier integral of 
a “conhning” gluon propagator 1 jk^ [12]. The spectrum ob¬ 
tained in [33] is in general agreement with the expected val¬ 
ues. As discussed here, this procedure is not able to generate 
a linearly rising potential, associated with conhnement in the 
static case (see also [19]). In our study, we have ensured that 
the expression in Eq. (21) approaches 1 jk^ in the ultraviolet 
limit by rescaling the corresponding term in the potential by 
a suitable constant. 

We note again that our aim was to gain a qualitative un¬ 
derstanding of the interplay between perturbative and non- 
perturbative features of the interquark potential. As veri- 
hed in our study, even though the full nonperturbative gluon 
propagator was used, the potential is non-confining, i.e. con¬ 
hnement is washed away by the use of the (tree-level) pertur¬ 
bative approximation for the interaction. Nevertheless, the 
resulting potential (with the addition of a linear term) pro¬ 
vides a slightly better description of the spectra, with the 
same number of ht parameters as the Cornell potential. 

The authors thank B. Blossier and F. Navarra for use¬ 
ful comments. W.S. thanks the Brazilian funding agencies 
FAPESP and CNPq for hnancial support. A.C. and T.M. 
thank CNPq for partial support. 
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Table 1 Experimental spectrum of charmonium and input values considered in our fits (see text). 


Particle Name 

Mass (GeV) 


/ 

Mass Input l(GeV) 

Mass Input 2(GeV) 

T7,{15) 

2.9836(7) 

3.096916(11) 

0-+ 

1— 

0 

3.06857(17) 

3.040(57) 

xcom 

xcim 

hc(lP) 

xcim 

3.41475(31) 

3.51066(7) 

3.52538(11) 

3.55620(9) 

0++ 

1++ 

1+- 

2++ 

1 

3.52528(8) 

3.485(70) 

77,(25) 

V/{25) 

3.6394(13) 

3.686109(14) 

0-+ 

1 — 

0 

3.67442(32) 

3.663(23) 

1/7(3770) 

3.77315(33) 

1 — 

0 or 2 

3.77315(33) 

3.77315(33) 

X(3872) 

3.87169(17) 

1++ 

1 

3.87169(17) 

3.87169(17) 

Xco(2P) 

Xc2{2P) 

3.9184(19) 

3.9272(26) 

0++ 

2++ 

1 

3.9257(25) 

3.9228(44) 

1/7(4040) 

4.039(1) 

1 — 

0 or 2 

4.039(1) 

4.039(1) 

1/7(4160) 

4.191(5) 

1 — 

0 or 2 

4.191(5) 

4.191(5) 

X(4260) 

4.251(9) 

1 — 

0 or 2 

4.251(9) 

4.251(9) 

X(4360) 

4.361(13) 

1 — 

0 or 2 

4.361(13) 

4.361(13) 

1/7(4415) 

4.421(4) 

1 — 

0 or 2 

4.421(4) 

4.421(4) 

X(4660) 

4.664(12) 

1 — 

0 or 2 

4.664(12) 

4.664(12) 


Table 2 Experimental spectrum of bottomonium and input values considered in our fits (see text). We include the unconfirmed state ?)(,{15) 
we need it to improve our results, as discussed in the text. 


Particle Name 

Mass (GeV) 

y^c 

1 

Mass Input l(GeV) Mass Input 2(GeV) 

J?i(lS) 

r(i5) 

9.3980(32) 

9.46030(26) 

0-+ 

1 — 

0 

9,4447(10) 

9.429(31) 

xmm 

Xhim 

hb{lP) 

Xhii^P) 

9.85944(42) 

9.89278(26) 

9.8993(10) 

9.91221(26) 

0++ 

1++ 

1+- 

2++ 

1 

9.89970(44) 

9.886(26) 

r(25) 

10.02326(31) 

1 — 

0 

10.02326(31) 

10.02326(31) 

r(iD) 

10.1637(14) 

2— 

2 

10.1637(14) 

10.1637(14) 

Xho{2P) 

10.2325(4) 

0++ 




Xbi{2P) 

10.25546(22) 

1++ 

1 

10.26022(43) 

10.251(18) 

Xb2i2P) 

10.268650(22) 

2++ 




r(35) 

10.3552(2) 

1 — 

0 

10.3552(2) 

10.3552(2) 

Xb{iP) 

10.534(9) 

??+ 

1 

10.534(9) 

10.534(9) 

r(45) 

10.5794(12) 

1 — 

0 

10.5794(12) 

10.5794(12) 

r( 10860) 

10.876(11) 

1 — 

0 or 2 

10.876(11) 

10.876(11) 

r(ii020) 

11.019(8) 

1 — 

0 or 2 

11.019(8) 

11.019(8) 



